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We study the spatiotemporal properties of coherent states (peaks, holes, and fronts) in a bistable activator- 
inhibitor system that exhibits biochemical saturated autocatalysis, and in which fronts do not preserve spatial 
parity symmetry. Using the Gierer-Meinhardt prototype model, we find the conditions in which two distinct pin- 
ning regions are formed. The first pinning type is known in the context of variational systems while the second is 
structurally different due to the presence of a heteroclinic bifurcation between two uniform states. The bifurca- 
tion also separates the parameter regions of counterpropagating fronts, leading in turn to the growth or contrac- 
tion of activator domains. These phenomena expand the range of pattern formation theory and its biomedical 
applications: activator domain retraction suggests potential therapeutic strategies for patterned pathologies, such 
as cardiovascular calcification. 

PACS numbers: 89.75.Kd, 47.20.Ky, 47.54.-r, 87.18.-h 



Reaction-diffusion (RD) systems display many kinds of 
pattern formation, ranging from diverse stationary patterns 
to breathing and/or replicating spots and spiral waves QJJ] . 
One important application of RD systems is to biological pat- 
tern formation and morphogenesis, as originally suggested 
by Turing: reaction and diffusion of two chemical sub- 
stances ("morphogens"), at different diffusion rates, can in- 
duce the spontaneous symmetry breaking of a homogeneous 
state and thus form patterns |2]. Turing described a sim- 
ple activator-inhibitor system; since then, many model equa- 
tions in physical, chemical, and biological contexts have been 
shown to form periodic patterns via a Turing-type instabil- 
ity |Q], S St . There is a significant difference between chemi- 
cal systems and biological media: In biology, the formation of 
labyrinthine patterns, for example, has been attributed to the 
phenomenon of saturation yi|5|]. 

Coherent states are another important property of spatially 
extended media 16J]; we distinguish between heteroclinic or- 
bits (fronts) and homoclinic orbits (localized states) in one- 
dimensional (ID) physical space. In the past, the formation of 
localized states in bistable systems was often linked to interac- 
tions between two front solutions |0] • However, this approach 
requires a large separation between two interacting fronts and 
thus is not valid for the analysis of spatial homoclinic orbits. 
By a different method, it was recently shown that localized 
states in variational systems are intimately related to Maxwell 
points, and thus nonpropagating fronts, connecting either two 
uniform states [8] or a uniform and a periodic ]9[] state. In 
the second case, there can be an effective broadening of the 
Maxwell point so that stationary fronts exist over a finite pa- 
rameter range 11011 : this phenomenon is referred to as pinning. 

While biological patterns have been extensively studied in 
the context of periodic patterns [3], there is no unified the- 
ory of coherent states (localized and fronts) and the resulting 
spatiotemporal dynamics in systems that involve autocatalytic 
saturation [11]. On the other hand, the effect of front asym- 
metry in the context of formation mechanisms of localized 
states is also unclear [6J] . To advance these issues, we study the 
formation and the interaction mechanisms of coherent states, 



holes, peaks and fronts, in a biological activator-inhibitor pro- 
totype model. We show that the known phenomenon of con- 
tinuous branch pinning (CBP) can also be found in dissipative 
systems. In addition, we present a type of pinning, to which 
we refer as discontinuous branch pinning (DBP). In case of 
CBP, the existence region of holes, which is below the Turing 
onset, is disjunct from the region of peaks, which is near the 
bistability saddle-node of the uniform states. As the control 
parameter is varied, the Turing onset approaches the saddle- 
node, and the regions of holes and peaks overlap. This over- 
lap results in the formation of a DBP structure of holes via 
the time-independent front (hetoroclinic bifurcation), which 
in turn also separates regions of counterpropagating asymmet- 
ric fronts. As a consequence, front counter-propagation and 
the existence of localized states give rise to an intriguing be- 
havior: activator domain growth or contraction accompanied 
by a simultaneous formation of embedded localized axisym- 
metric or striped states, as shown in Fig.[TJ 

We start with a phenomenological activator-inhibitor model 
introduced by Gierer and Meinhardt IU2I1 . to study biological 
pattern formation. Our interest in this model stems from calci- 
fication in vascular-derived mesenchymal stem cell (VMSC) 
cultures lfl3ll . cells that are considered to be responsible for 
cardiovascular calcification 11411 . In particular, it is assumed 
from biochemistry that the activator obeys a saturated au- 
tocatalytic reaction lfl5ll and promotes production of the in- 
hibitor [16]. The cells also express a rapidly diffusing in- 
hibitor [13]. This biochemical activator-inhibitor dynamics 
therefore qualitatively rationalize the framework of Gierer- 
Meinhardt model d, [lal . which in its dimensionless form 
reads fl7ll 
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where u(x,y) is activator, v(x,y) is inhibitor, D is the dif- 
fusion ratio, P represents the generalized cross reaction rates 



2 




t = 60 



r = 680 



t = 4000 




FIG. 1: (Color online) Top panel: Bistability region, S < Ssn, of 
uniform solutions (0, Wo) and (u+, v+) in the E — S plane (shaded 
region). The dark shaded region, St < S < Ssn, marks the lin- 
ear instability of the uniform (u+ , v+ ) state to nonuniform perturba- 
tions, where St is the onset of Turing instability. The dashed line 
S = Sm, for E > E c ~ 3.54, represents the stationary front 
between (0, vo) and (u+,v+) states (heteroclinic bifurcation). For 
S < Sm, («+,«+) invades (0, vo) while for Sm < S < St (0, vo) 
invades («+, Bottom panel: Demonstration of the activator (a) 
expansion and (b) contraction process at marked locations in the left 
panel, stationary front S = Sm — 3.35 for E = 9, via numeri- 
cal integration of Eq. (TJ on x = y = [0, 15] domain with periodic 
boundary conditions. The frames show a gray-scale map of the it 
field, where black denotes a high value of u. The initial condition is 
either uniform (a) u = 0.01 or (b) u = 1 in the left half of the do- 
main, while the right half is a uniform random distribution between 
0.01 and 1. Parameters: (a) S = 3.3, (b) S = 3.4, and for both 
P = 1, D = 0.005. 



ratio, E is the degradation ratio and S is the generalized in- 
hibitor source term. Due to the interest in vascular calcifi- 
cation [14], we chose here the parameters estimated in lfl3ll 
(see supporting online supplement), i.e., P ~ 0(1), D ~ 
0(1O -2 — 10 -3 ), and focus on the two remaining parameters 
S and E, which are experimentally accessible. We note that 
localized peaks have also been studied in other forms of Eq. 
®, in ID and 2D HUH. 

Equation (Q]i admits three uniform solutions, one of which 
(u, v) — (0, vo) = (Q,S/E) is trivial, and two nontrivial 
(u,v) = (u±,v±) IU9I1 . Uniform states appear as 

unstable states from S = and annihilate with the stable 
(u+,v + ) branch at a saddle-node bifurcation, S = Ssn, cre- 



ating asymmetric bistability between (0, vq) and (u + , v+) for 
< S < Ssn- A primary requirement for both stable local- 
ized and front solutions is the linear stability of the uniform 
solutions to nonuniform perturbation. Standard linear analy- 
sis of (0, vo) and (u+, v+) yields that only the nontrivial state 
goes through a Turing instability at S = St- The unstable 
Turing region St < S < Ssn shrinks by approaching the 
saddle-node bifurcation, S = Ssn, as the degradation ratio 
E is increased (see Fig. [TJ. Consequently, stable front solu- 
tions connecting the uniform (u+,v+) and (0, Vq) states exist 
below S = St- 

Next, we study the properties of localized states via spa- 
tial dynamics 12011 followed by numerical branch continuation 
method 1I2TI1 and temporal eigenvalue numerical analysis to 
determine stability. We set dtu — dtv = and analyze ([TJ) 
as a reversible fourth order ordinary differential equation due 
to the x — > —x invariance. Linearizations about the uniform 
states result in steady state solutions which are proportional to 
exp(Aa;), where the spatial eigenvalue A satisfies a fourth or- 
der algebraic equation. For the (u + , v+) case, three possible 
solutions for A in complex phase space arise: (J) for S > St 
the eigenvalues split on the imaginary axis; (if) for S < St 
they split and form a complex quartet; and (fff) at the Turing 
onset, S = St, a double multiplicity of pure imaginary eigen- 
values A = ±ikT is present, where fcx is the critical Turing 
mode. Thus, localized states that exponentially asymptote to 



can be expected to arise for S < Sj- (9|. On the 



other hand, linearization about (0, vq) 



yields four real eigen- 
Thus, localized states 



values, A = ±\/E and A 
that exponentially asymptote to (0, «o) are expected to arise 
everywhere. For finite values of E, effective bifurcation of 
these states is at the transcritical onset, 5 = (the mechanism 
will be discussed elsewhere). In the following, we refer to L + 
as "holes" while to L° as "peaks" in the background of the 
nontrivial and trivial activator fields, respectively. 

To study the localized solutions far from the onset, we im- 
plement a numerical branch continuation method [21]; the re- 
sults are presented as a function of S in terms of the norm 
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where L is the spatial period; the domain size used in our 
numerical calculations is much larger than the spatial period 
2-k/Ut- We distinguish between two qualitatively distinct re- 
gions (see Fig.[TJ: E < E c , where the holes and peaks do not 
overlap, and E > E c (see Pig. EJ, where peaks and holes re- 
gions do overlap, and reconnect at S — Sm via the stationary 
heteroclinic cycle between v+) and (0, Vq) (see Fig. [3). 

First, we discuss the E < E c case. Localized hole solu- 
tions, L + , and periodic Turing states (not shown) bifurcate 
subcritically as unstable small amplitude states 11711 . from the 
Turing onset, S — St, and form a CBP region (shaded re- 
gion in Fig. 0, as also happens in variational systems [9]. 
Since the localized holes bifurcate from a nontrivial state, only 
two even parity solutions with phase shifts of n can form. 
Thus, we termed a branch with an odd number of holes as Lq 
[Fig. |2 a)] and the one with an even number as L+ [Fig.[2fb)]. 
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FIG. 2: (Color online) To 
ID, showing the norm © 
and localized states: Lq , 
stable (unstable) solution 
region. St and Ssn cor 
(see text for details), resj 
(dashed) profiles at locati 
The saddle-nodes and thi 
branches agree with the st 
(a,b) S ~ 0.071 (c) S ~ 
D = 0.005, P = 1. 
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Each stable branch of odd and even hole solutions (inside 
the shaded region) indicate increasing (with N) number of 
holes [9]. Here we used a relatively small domain which lim- 
its further excursions of the two branches (Fig. |2j, while on 
an infinite domain, the branches contain an infinite number 
of saddle-node points. On the other hand, the localized states 
that are homoclinic to the (0,vq) state effectively bifurcate 
from S — and do not form a pinning region, due to real 
eigenvalues. These states contain only a single stable branch 
of Lq while the other even branch L\ is unstable [see profiles 
(c,d) in Fig. |3, denoting a repulsive interaction between two 
neighboring peaks. Here the subscripts and ir, are artificial 
notations. 

Next, we turn to a discussion of the second case, E > E c . 
We showed in Fig. [TJ that an increase in E results in St — > 
Ssn- Above E = E c ~ 3.54, this shift causes an overlapping 
between the CBP region of hole states and the peak states that 
are present near the saddle-node S — Ssn- This hole-peak 
interaction allows formation of a time-independent asymmet- 
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FIG. 3: (Color online) Top 
panel: Bifurcation diagram 
for E = 4 (solid (dashed) 
lines mark stable (unsta- 
ble) solutions). The inset 
in the top panel clarifies 
the behavior of Lq and Lq 
branches at S = Sm — 
1.233. Bottom panels: Pro- 
files of u at locations in- 
dicated in the inset. Note 
the asymmetry of the het- 
eroclinic orbit (half of the 
cycle) in (c). Other param- 
eters as in Fig. [2] 



ric front (via heteroclinic bifurcation) at S = Sm, connecting 
the trivial and the nontrivial states (see Fig.O. All the local- 
ized state branches that were present in the E < E c case, now 
effectively collapse nonmonotonically at S = Sm- The single 
and the double hole state, L = L + , still bifurcate from the 
Turing onset (see Fig. [3]), but all other branches (not shown 
here for simplicity) arise from and reconnect after a single ex- 
cursion, back to S = Sm- This DBP region is structurally 
different from the classical CBP, described for E < E c and 
in 09J. The L° branches also go through a significant change: 
while L° remains unstable, the Lq branch exhibits decreasing 
oscillations around S = Sm (see inset of Fig. [3]). This behav- 
ior is known to give rise to additional localized states [8, 22], 
with weak spatial oscillations around the primary peak [see 
Fig. Eft)]. 

The activator expansion or contraction and formation of 
isolated stripes and axisymmetric localized states (Fig. [TJ is 
now transparent. The effective merging point, S = Sm, sep- 
arates regions of counterpropagating single fronts. Thus, op- 
posite front velocities and the simultaneous presence of stable 
localized holes, peaks and stripes, result in a nontrivial behav- 
ior for bistable nonvariational systems with a broken parity 
symmetry of fronts [see Fig. [3jc)]. In the absence of local- 
ized states, the system converges to a finite uniform activator 
state, a phenomenon that is qualitatively similar to coarsen- 
ing in variational systems [7]. An additional important aspect 
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is the stability of coherent structures to curvature and trans- 
verse perturbations in 2D. While stability thresholds of local- 
ized states may indeed change in 2D, the ID analysis delimits 
the regions where stable localized states can also be obtained 
in 2D ifni l22tl . By contrast, instability to transverse pertur- 
bations of the front solutions or large circular domains may 
lead to the formation of periodic patterns that resemble peri- 
odic patterns arising from the Turing instability [23]. In our 
parameter range these fronts are transversally stable, but we 
cannot exclude the existence of transverse instabilities in gen- 
eral. 

We have demonstrated a mechanism for the expansion and 
contraction of activator domains, and the respective formation 
of localized states, in a bistable activator-inhibitor system 
with saturated autocatalysis (Fig. [TJ. This mechanism stems 
from the formation of an asymmetric heteroclinic orbit (at a 
critical parameter value), merging the branches of both hole 
and peak states. Turing instabilities are found to be important 



for such behavior, and also for a discontinuous branch 
pinning structure, as shown in Fig. [3] Since we used global 
bifurcation methods, we expect similar scenarios to arise in 
other bistable systems as well, precisely for the same reason 
that pinning [9] is found in variational models and in our 
dissipative system. Furthermore, our theoretical results sug- 
gest a distinct spatiotemporal behavior in biological systems 
exhibiting saturation [5]. The concept of activator contraction 
has potential biomedical applications, for example, localized 
activator-induced "nodules" are common in biomedicine, 
such as spotty formations in atherosclerotic calcification 112411 . 
The ability to produce conditions of "activator contraction" 
to reduce or eliminate these localized spots, therefore has 
intriguing therapeutic possibilities. 

We thank L.L. Demer, Y. Tintut, K. Bostrom, E. Knobloch, 
and J. Burke for stimulating discussions. 
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